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1 Liebmann technical documentation series 


1. Wyznaczanie rozktadu pola elektrostatycznego w prézni metodg relak- 
sacyjna Liebmanna. (Polish version / wersja polska) 


2. Determination of electrostatic field distribution by using Liebmann relax- 
ation method. (English version / wersja angielska) 


3. Graphics. Mapping voltages to colours. (colormaps) 


4. Laplace equation 2D (XY). (Cartesian coordinates). Relaxation scheme 
explained. (5 - point star) 


5. Laplace equation 2D (ZR) (Cylindrical coordinates). Relaxation scheme 
explained. (5 - point star) 


6. Liebmann source sode. (ANSI C programming language) 


2 Versions of this document 
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. version 3 - 2024.02.02 
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. version 5 - 2024.05.18 
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. version 7 - 2024.05.24 


. version 8 - 2024.06.06 (complete P;..Po) 
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. version 9 - 2024.06.09 


3 Solving Laplace equation using relaxation method 


| tried to solve Laplace equation using mainly information from Pierre Grivet’s 
book (Electron Optics) - [1]. 
There are few editions of this book (1965, 1972). Second edition (1972) con- 
tains explanation of relaxation method (page 38). 

More generalized approaches has been drafted by James R. Nagel - [2]. 


https://my.ece.utah.edu/~ ece6340/LECTURES/Feb1/ (visited 2023-03-01). 
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Taylor expansion in cylindrical coordinates has been found on the Internet: 
3}. 


There are also publications edited by Albert Septier: Focusing of Charged 
Particles [4] and Applied Charged Particle Optics (part A). [5]. 

| have also found some ideas in publication of D W O Heddle: Electrostatic 
Lens Systems [6] (especially using PC computers to solve electrostatic prob- 
lems). 

| have also found (brief) description of by - hand solving of Laplace equa- 
tion by Bohdan Paszkowski - [7] (Polish edition). English translation of this book 
also exists - [8]. 


| would like to thank many people, who helped me with this challenge. Espe- 
cially prof. dr hab. Mieczystaw Jatochowski (Supervisor of my master’s thesis), 
who enabled me to use SIMION and MATLAB software while writing master’s 
thesis about electron optical systems at University of Maria Curie - Sktodowska 
in Lublin in 2008. | would also thank to prof. Marcin Turek for fruitful discus- 
sion about numerical methods. What is more, my colleague Bartosz in 2012 
had explained me general problems with software efficiency. So he had also 
contributed significantly to the idea of Liebmann software (especially using C 
language). 


is 4 Explanation of symbols in calculations 


169 ¢ P, - i-th mesh node 

170 Vj - value of electrostatic potential at node P;. Unit - [V] 

171 * h-mesh step (for example h,, - mesh step in x direction). Unit - [mm] 
172 * gi+/—- gradient in direction i (for example gi-— = Yao Fiee Unit - [|] 
173 * drow - index of row in mesh. Values of 7,5, = 1, 2, .., size_row 

174 * ico - index of column in mesh. Values of 7..; = 1, 2, .., size_col 

175 * p-in book: - [1] r = ph,., so for off - axis point we have: p = (tow — 1) 
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5 Mesh ZR - type A (on axis) 


hig Fhe 
gradient V outside a mesh exists 
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Figure 1: Mesh ZR type A 


179 


180 


181 


6 Mesh ZR - type B (on axis) 
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Figure 2: Mesh ZR type B 


2 7 Mesh ZR - type C (on axis) 
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Figure 3: Mesh ZR type C 
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8 Mesh ZR - type D (on axis) 
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Figure 4: Mesh ZR type D 
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9 Example of A-type mesh in ANSI C (on axis) 


Example of A- type mesh in ANSI C program. The mesh is represented by 2 
dimensional array of double precision numbers. Rows and columns in mesh 
are numbered from 1 (this was my choice) instead of default 0 (as usual in C 
language). This choice nas pros and cons. Is is easier to calculate mesh size 
(size_row * size_col). Access to each node can be also more intuitive, but logic 
in each library function must contain this shift between node ordering styles. 


—+| 0 1 2 3 4 5 
Gr+ 
trow 
4 3 18 19 20 21 22 23 3 
3 2 12 13 14 15 16 17 2 
2 1 6 i 8 9 10 11 1 
1 0 / 0 1 2 3 4 8) 0 
4 | 
Gz Get 
—| 0 1 2 3 4 5 
Gr— 


1 2 3 4 5 6 tol 


Figure 5: ANSI C - mesh XY type A 


Note. This is more general example of off-axis” mesh. If bottom egde of 
mesh lies on axis Oz, then gradient g,— does not exist. 


* gz— =double* ptr_gZ_minus 
* G24 =double* ptr_gZ_plus 
* gr— = double* ptr_gR_minus 


* Or+ =double* ptr_gR_plus 
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¢ V = double* ptr_V 


unsigned int size_row == 4 


unsigned int size_col == 6 


unsigned int i_row == 1, 2, .., 4 


unsigned int i_col == 1,2, .., 6 


double h_z == 1.0 [mn] 


* double h_r == 2.0 [mn] 


The following picture describes analogous version of ptr_V mesh, which 
can be dynamically allocated on heap by pointer metod. The mesh is rep- 
resented by single block of memory. The numbers or rows and columns are 
also known, so each node can be also accessed by appropriate index (memory 
address). 

The following picture describes analogous version of mesh, which can be 
easily dynamically allocated on heap by pointer metod. The mesh is repre- 
sented by single block of memory. The numbers or rows and columns are also 
known, so each node can be also accessed by appropriate index (memory ad- 
dress). 


ptr_v—— | 0 1 2 Ag 23 


Figure 6: ANSI C - mesh ZR type D 


Each mesh point has its unique index (let’s say icp - (index of central 
point)), which can be determined, if we know indices of row and column (i_row, 
i_col). 

icp == ( i_row - 1) * size_col + i_col - 1 (9.1) 


For example for each point of a mesh indices of row and column have val- 
ues: 


> $ize_row 


i] 
ll 
rae 
. 
. 


1i_row 
(9.2) 


i_col > $ize_col 


i] 
ll 
rae 
. 
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223 Also for any relaxation formula for off - axis case the p symbol appears. This 
224 symbol is connected with r cylindrical coordinate of given node: 


225 


r= ph, (9.3) 


226 SO: 


p == ( i_row - 1) (9.4) 
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10 Example of B-type mesh in ANSI C (on axis) 


Example of B- tyoe mesh in ANSI C program. The mesh is analogous to A - 


type mesh. There are no electric field gradients on mesh borders. 


trow 
4 18 19 20 21 22 23 
3 12 13 14 15 16 17 
2 6 7 8 9 10 11 
1 / 0 1 2 3 4 i) 
V 
1 2 3 4 i) 6 


Yeol 


Figure 7: ANSI C - mesh XY type B 


V = double* ptr_V 
unsigned int size_r 
unsigned int size_c 
unsigned int i_row 
unsigned int i_col 
double h_z == 


double h_r == 


Ow 


ol 


1.0 [mm] 


2.0 [mm] 


== 4 
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Example of C- type mesh in ANSI C program. The mesh is analogous to A - 


Example of C-type mesh in ANSI C (on axis) 


type mesh. Just mesh mesh step hz = h, = h. 


— | 0 1 2 3 4 8) 
Gr+ 
trow 
4 3 18 19 20 21 22 23 3 
3 2 12 13 14 15 16 17 2 
2 1 6 7 8 9 10 11 1 
1 0 j 0 iL 2 3 4 5 0 
I: 
Gz Get 
— | 0 1 2 3 4 8) 
Gr— 
1 2 3 4 5 6 teol 


Note. This is more general example of ,,off-axis” mesh. If bottom egde of 
mesh lies on axis Oz, then gradient g,_ does not exist. 


* gz— =double* ptr_gZ_minus 


Figure 8: ANSI C - mesh XY type C 


gz+ = double* ptr_gZ_plus 
gr— = double* ptr_gR_minus 


Gr+ = double* ptr_gR_plus 


V = double* ptr_V 


unsigned int size_row == 
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unsigned 
unsigned 
unsigned 


double h 


int size_col 
int i_row == 
int i_col == 


== 1.0 [im] 
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22 12 Example of D-type mesh in ANSI C (on axis) 


253 


Example of D- type mesh in ANSI C program. The mesh is analogous to B - 


asa type mesh. Just hy = hy =h. 
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Figure 9: ANSI C - mesh ZR type D 


V = double* ptr_V 


unsigned 
unsigned 
unsigned 
unsigned 


double h 


int size_row 
int size_col 
int i_row == 
int i_col == 


== 1.0 [im] 


== 4 
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x1 13 Partial derivatives on Oz axis 


22 13.1 Personal note 

263 This is my presonal interpretation. | cannot guarantee correctness of this ap- 
264 proach 

25 13.2 Nodes numbering (on axis O.) 


2e6 Wewill try to work with P2 point (determine approximations of aprtial derivatives 
267 for point P2, which lies on axis O,). Nodes numbering on axis O, differs from 
2e8 Numbering convention in Pierre Grivet’s book. 


Figure 10: Nodes on axis Oz 
269 Point P2 is situated on O, axis. It has 2 neighbours on axis O, : P, and P3. 


270 Node Ps lies abowe P, node. The mesh step in r direction is h,. The mesh 
271 Step in z direction is h,. 


2 13.3 Taylor expansion of V(.,,) function in cylindrical coordinates 


27a Taylor expansion of V/,,.. function in cylindrical coordinates (3): 


OV 
Ver) = Veco + ee (z _ z)+ 


20,70) 
OV 
Or ( ae (13.1) 
z0,70 
1 (=) ‘ 
~ | > (z— z0)° +... 
21 ax (20,70) 


274 13.4 Laplace operator in rotationally symmetrical systems 


275 Laplace operator in cylindrical coordinates (cylindrical symmetry) consist of 3 
27s elements [1] (on page 42): 
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OV 
V* (Yen) = (5) 7m 


: (>) + (13.2) 


r\ or 
OPV 
Oz? 
277 In this chapter we will try to determine approximation of each term. 


2 13.5 Value of first partial derivative of V with respect to r on axis 
279 Oz 


280 In cylindrically symmetrical field first partial derivative of V (by r) on axis Oz 
ze1 equals zero (because V(,4,) = Vi—ar)) 


(F) = (13.3) 
or (z,r=0) 


222 13.6 Value of second partial derivative of V with respect to , on 
283 axis Oz 


zea In this subchapter we will try to determine the first term of equatior/13.2] 

285 In our case there is node P, on axis Oz. The nearest neihgbour of P is 
28s node P;, which lies ,over Oz axis”. The distance between P: and P; is h,. 
227 When we ,, walk away” axis Oz in r direction (from point P, to point Ps), the 
2sa_ electric potential V; can be determined from truncated Taylor expansion |13.1 
29 by expression: 


1 2. 
V5 2 Vo+ (¥) hp + — (Sr) he? (13.4) 
Or Py 2! \ Or Py 
290 We want to determine the second derivative: 
2 
as =? (13.5) 
Or? Ps 


291 We solve equation (using relation 


5 2! (Vs — 2 (V5 — 
OV)  2(Vs— Va) _ 2(V5 — Va) (13.6) 
Or? } p, h2 Ae 
292 This is final form of approximation of the second derivative of V with respect 


23 to r on axis Oz. It will help us to determine Laplace operator in rotationally 
24 Symmetrical systems. 
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25 13.7 Value of first partial derivative of V with respect to , divided 
296 by r on axis Oz 


297 We will try to determine the second term of relation [13.2] Wheh we are on Oz 
298 axis, the second term has to be determined (because it aims to value lf ). 

299 When we ,, walk away” axis Oz in r direction, the electric potential V(., ) 
300 Can be determined from truncated Taylor expansion by: 


301 


OV 
Vizo,r) © V(zo,ro) + (+) (r — 19) (13.7) 
® J (20,70) 
302 On Oz axis ro = 0, SO (r9 — TT) =r 
303 
304 Thus we have: 
OV 
V XV + (+) -T (13.8) 
(20,7) (20,0) 
: ° Or (z0,0) 
305 Now let us differentiate (both sides) of such relation: 
O 
a 13. 
le (13.9) 
306 We get: 


2 
(>) ~() ‘(F) r+ (ZF) 1 (13.10) 
OF J tea) T } (20,0) Te T J (20,0) 


307 On axis Oz we can apply relation That’s why we can remove these 
30s two terms (first and third) from equation|13.10 
309 So we get (if 7 = 0): 


2 
(+) ~ (= _ T (13.11) 
Or (zo,r) Or (z0,0) 
310 We can now divide both sides by r. 
1 
| == (13.12) 
r 
311 We have relation, which has been published in Pierre Grivet’s book|1]. 
1 ed 
(=) “5 (= >) (13.13) 
r Or on Or (xa0) 
312 Approximation of this term on numerical mesh has been already determined 


sis in previous subsection (13.6). 
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13.8 Value of second partial derivative of V with respect to z on 
axis Oz 


The third term of Laplace operator in rotationally symmetrical systems |13.2 
takes form ( on picture|10}): 


2 Va—V2 _ V2-Vi 
—2 
(=) ~ he he = rs VW (13.14) 
Po 


hz h2 


Zz 


Now we have determined all the 3 approximations o partial derivatives of V 
in cylindrically symetrical systems (on axis O,). 
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14 Partial derivatives off Oz axis 


14.1. Personal note 


This is my presonal interpretation. | cannot guarantee correctness of this ap- 
proach 


14.2 Nodes numbering in Liebmann mesh (off axis O.) 


A 
Py Pg Po 
r O O O 

Py, Ps Ps 

O @ O 4 

hp 

P, P. P. 

a on a T5 = ph, 

hz 
Oz 


Figure 11: Nodes off axis Oz. Exemplary vector r5 describes distance from 
axis O, to node P; 


Mesh step in z direction is h,. Mesh step in r direction is h,. Sample mesh 
points Ps lies off O, axis. Distance between mesh point P; and O, axis is rs. 

For ANSI C meshes (in Liebmann source code) the following relations have 
place: 


r = ph, (14.1) 


(14.2) 


Dp =trow — 1 


14.3 Taylor expansion of V(.,.. function in cylindrical coordinates 


Taylor expansion of V, ,) function in cylindrical coordinates (3): 
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OV 
Vizr) => Vizo,ro) + ey) (z = zo)+ 
20,70 


ov 
Op (r —r9)+ (14.3) 
(20,r0) 


1 (S) ‘ 
~|a=> (z— zo)" +... 
ALN On? 7 (ear6) 


331 14.4 Laplace operator in rotationally symmetrical systems 


332 Laplace operator in cylindrical coordinates (cylindrical symmetry) consist of 3 
333 elements [1] (on page 42): 


OV 
V? (Very) = (Sr) 1 
1 (OV 
: (> ) i (14.4) 
(3) 
Oz? 


334 In this chapter we will try to determine approximation of each term. 


33 14.5 Value of second partial derivative of V with respect to 7 off 
336 axis Oz 


(14.5) 


= hy h2 
337 14.6 Value of first partial derivative of V with respect to , divided 
338 by r off axis Oz 


1 @ _1-e _ V—-Ve 
Ps 


— Re = 14. 
Or ry 2hp 2rshz ee) 


T5 


30 14.7 Value of second partial derivative of V with respect to z off 
340 axis Oz 


& _ Ve Ms _ V5 V4 7 Yt = 9 


~ (14.7) 
A hz h2 
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15 Relaxation formula for node P1 (on axis Oz) 
15.1 Node description 


Left, bottom corner of mesh ZR (on axis Oz). 


15.2 Calculation of relaxation formula 


Laplace equation at node P; 


OV, 10V; OV, 
(z,r) Eay (z,r) + (z,r) =§ 
Or? Et Or Ps Oz? 
Py 1 Pi 


Approximation of partial derivatives of Vi,,-) at node P;, 


(32) ae ae Vi) 
Py 


Or? - hy RB 


10V(z.7) nf, 2 (V4 _ Vi) 
r Or Jp h2 


V: —V, Vi -V; 2 
Ver) tay a he _-VN ge 

Oz? = hz h2 hz 
Let us substitute approximations to Laplace equation. 


2(Va—Vi) , 20Va—-Vi) Va-M _ me- 


he ho ORR a 


2(V4—Vi) | 2(0V4a—-V) | Vo-Mi _ 91z- 


h2 h2 / Re hz 


Let us find Vj 


Y=? 


Let us multiply both sides 
| -AZh; 


We obtain 


QVih2 — 2V,h2 + 2Vyh2 — 2h? + Veoh? — Vib? = 
Qiz—hzh? 
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(15.1) 


(15.2) 


(15.3) 


(15.4) 


(15.5) 


(15.6) 


(15.7) 


(15.8) 


(15.9) 


(15.10) 


352 


353 


354 


355 


356 


357 


Let us simplify this equation: 


V, (2h2 + 2h2 +h?) = 
Qh? + 24h? + Voh2 — giz—hzh? 


So we have: 


V, (4h2 + h2) = AV gh? + Voh2 — giz—hzh? 


15.3 Final forms of relaxation formula 
15.3.1 zrLV_RELAX5 P1_ON A 
hz # hy 
Nz- FO 


ee AVgh2 + Voh? — giz—hzh? 
oe 4h + h2 


15.3.2 zrLV_RELAX5 P1_ON B 
he # hy 
Qz- = 0 


ee AVgh? + Voh? 
1 Gh? + h2 


15.3.3 zrLV_RELAX5 P1_ON_C 
hz=hp=h 
Jiz— = 0 


— 4V4 + V2 — giz-h 
ei 5 


Vi 


15.3.4 zrLV_RELAX5 P1 ON D 
hz, =h,=h 
Nz = 0 


— 4V4+ V2 


Yi 5 
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(15.11) 


(15.12) 


(15.13) 


(15.14) 


(15.15) 


(15.16) 


ss 16 Relaxation formula for node P2 (on axis Oz) 
30 16.1 Node description 


3c0 Bottom edge of mesh ZR (on axis Oz). 


31 16.2 Calculation of relaxation formula 


3e2 Laplace equation at node P» 


2 
V (View) p, =0 (16.1) 
PV. r) i ae r) PV, r) 
( D2 + (C2 ). + 722 =0 (16.2) 
P»2 2 P2 
363 Approximation of partial derivatives of V(,,) at node P, 
Ves) > tots -_ a 2 2 (V5 = V2) He _ 
Or? - hy RP 
P2 
1LOViGs 2 (V5 — 
COR, Pe he 
O7View e Laat _ wo ” VY, + V3 — 2V2 (16 5) 
Oz hs hg 
P2 
364 Let us substitute approximations to Laplace equation. 
2(V5— V2) | 2(Vs—Va) | Vit V3— 2V2 
72 2 2 =0 (16.6) 
365 There are no g expressions to move, to formula 7 has identical form as 
see formula 6. 2(Vs— V9) 2(Vp—Va) . Vi + Va — 2V 
eee ee) OS aes eee ee 
72 + 72 12 =0 (16.7) 
367 Let us find Vo 
Y=? (16.8) 
368 Let us multiply both sides 
| -Azhe (16.9) 
369 We obtain 


QV5h2 — Woah? + 25h? — Woh? + Vih?2 + Vzh?2 —2Voh2 =O (16.10) 


26 


370 


372 


373 


374 


375 


376 


Let us simplify this equation: 


Vo (2h2 + 2h? + 2h?) = 2V5h?2 + 2Veh2 + Vib? + Vgh? 


So we have: 


Vo (4h? + 2h2) = 4V5h2 + (Vi + V3)h2 


16.3 Final forms of relaxation formula 
16.3.1. zrLV_RELAX5 P2 ON A 
hz & hy 


_ AVaRZ + (Vi + Va)he 


V2 Ah2 + 2h2 


16.3.2 zrLV_RELAX5 P2 ON B 
hz & hy 


_ AVsh?2 + (Vi + V3)h2 


Va Ah2 + 2h2 
16.3.3 zrLV_RELAX5_P2_ON_C 
hehe =f 


— 4V5 + Vi + V3 


Vo 6 


16.3.4 zrLV_RELAX5 P2 ON D 
hz =hp =h 


_ 4V5+M4+ V3 


Vo 6 


27 


(16.11) 


(16.12) 


(16.13) 


(16.14) 


(16.15) 


(16.16) 


377 


378 


379 


380 


381 


382 


383 


384 


385 


386 


17 Relaxation formula for node P3 (on axis Oz) 
17.1 Node description 


Right, bottom corner of mesh ZR (on axis Oz). 


17.2 Calculation of relaxation formula 


Laplace equation at node P; 


ON GO: af SONG Ya hO Vien \, 23% 
Or? fOr Ps Oz? 
P3 7 P3 


Approximation of partial derivatives of V(,,) at node P3 


(52) stag est 2) 
P3 


Or? = hy RB 


(+e) a 2 (Ve = V3) 
P: 


r Or h2 
V32 —V: 3 
PV2n) 7 3 = 3 Va-Va _ Vo — V3 _ IBz+ 
Oz2 hz hz oh, 
P3 
Let us substitute approximations to Laplace equation. 
2(V6 — V3) 2(Ve—Vs) | Va- Vs , 93c+ _ 9 
h2 h2 © Rh, 
2(V6 — V3) | 2(Ve— V3) | Va— V3 _ get 
h2 h2 © Re hz 
Let us find V3 
Vas! 


Let us multiply both sides 
| -AZh; 


We obtain 


2Veh2 — 2h? + Wehz — Wahz + Vohe — Vahz = 
~9324h,h2 


28 


(17.1) 


(17.2) 


(17.3) 


(17.4) 


(17.5) 


(17.6) 


(77) 


(17.8) 


(17.9) 


(17.10) 


387 


388 


389 


390 


392 


393 


394 


Let us simplify this equation: 


V3 (2h2 + 2h2 +h?) = 
QWeh2 + 2eh2 + Veh? + g3z4hzh? 


So we have: 


V3 (4h2 + hz) = 4Veh2 + Voh? + giz—hzh? 


17.3 Final forms of relaxation formula 
17.3.1. zrLV_RELAX5_P3_ON_A 
hz # hy 
932+ FO 


Pe AVGh? + Voh?2 + g324hzh? 
al 4h2 + h2 


17.3.2 zrLV_RELAX5 P3 ON B 
igh 
93z+ = 0 


y, — 4Veh2 + Vohi 
3 AR? + he? 


17.3.3 zrLV_RELAX5 P3_ON_C 
he=hp=h 
J3z+ = 0 


— AVe + V2 + g3z+h 


ps 5 


17.3.4 zrLV_RELAX5 P3_ON _D 
h,=h,p=h 
932+ =9 


4Ve6 + V: 
a oe 


29 


(17.11) 


(17.12) 


(17.13) 


(17.14) 


(17.15) 


(17.16) 


395 


396 


397 


398 


399 


400 


401 


402 


403 


404 


405 


406 


18 Relaxation formula for node P4 
18.1 Node description 


Left edge of mesh ZR. 


18.2 Calculation of relaxation formula 


Laplace equation at node P, 


PV. r) iL ae r) PV, r) 
> + _ 5 ’ + 5 5 = 0 
Or 2 r Or P, Zz ce 
Approximation of partial derivatives of V(,,,) at node P, 


(note: r = ph, (such as in Pierre Grivet’s book) ) 


OV ese) 2 VWeaV4 — Va-v ee A ee eee 1 
OF? vs ~ hy re h2 
4 


fOr pf Die = 2ph? 


(=) oe ea V7- Vi 
P4 


= a. 
PV) 4. —— hae _Ve—V4 — 9az- 
Oz2 P h. h2 hz 

4 


Let us substitute approximations to Laplace equation. 


VitVr-2V4  Ve-VN . Vs—Va — gaz- 


he "Oph2 he hg TP 
VitVr—-2V4  Ve-VM . Vs—Va _ 9az- 
h2 ' Qph? heh, 
Let us find V4 
V4 =? 


Let us multiply both sides 


We obtain 


30 


(18.1) 


(18.2) 


(18.3) 


(18.4) 


(18.5) 


(18.6) 


(18.7) 


(18.8) 


(18.9) 


407 


408 


409 


410 


412 


413 


414 


415 


416 


QpVi hz + 2Vrh2 — ApVah2 + VzhZ — Vihz+ 
+2pVeh? — 2pVih? = 2pgaz—hzh? 


Let us simplify this equation: 


V4 (4ph + 2ph2) = Vi (2ph2 — h2) + V7(2ph2 + h2) + Vs2ph2— 


2pgaz—hz he 


So we have: 


Va2p(2h2 +h?) = Vih2(2p—1) + Vph2(2p+ 1) +Vs2phy — 2pgaz—hzh> 


18.3 Final forms of relaxation formula 


18.3.1 zrLV_RELAX5 P4_A 


hz & hy 
J4z— # 0 
h2(2p — 1) hz (2p + 1) Le g4z—hzhp. 
Va 2 a» Vit 2 a» 74 2 3 M5 2 2 
2p(2hz F h2) 2p(2hz ae h2) 2hz a he 2hz a hy. 
18.3.2 zrLV_RELAX5 P4 B 
he # bp 
gaz— = 0 
_ hQp-1) ),  beQptl) yy, 
4 Op(2h2 +2) 1" Qp(2h2 +2) '" 2h2+h2 ” 
18.3.3 zrLV_RELAX5 P4 C 
ite 
J4z— Fx 0 
2p —1 2p +1 1 gaz—h 
Mie Vi 4 Vt 
o> Bye er apg a 
18.3.4 zrLV_RELAX5 P4 D 
hy —hp=h 
Gaz— = 0 
2p —1 Mp +1 1 
Ve ey ee ey 
6p 6p 


(18.10) 


(18.11) 


(18.12) 


(18.13) 


(18.14) 


(18.15) 


(18.16) 


417 


418 


419 


420 


421 


422 


423 


424 


425 


426 


427 


428 


429 


19 Relaxation formula for node P5 
19.1 Node description 


Inner node of mesh ZR. 


19.2 Calculation of relaxation formula 


Laplace equation at node P; 


OV ies 10Viz7 PV r 
(z,r) Py pera d 4 (zr) _9 
Or? Or yp Oz? 
Ps > Ps 


Approximation of partial derivatives of V(,,,) at node Ps 


(note: r = ph, (such as in Pierre Grivet’s book) ) 


Tr 


hy h?2 


e 


(32) x es ee _ Vat Va — 2V5 
Or? ad re 
Ps 


1 V6.9) 2, UVeeeVs Vg — V2 
ae es 2ph? 


r Or ~ 


Oz? hz h2 


z 
Let us substitute approximations to Laplace equation. 


Vot+Ve—2V5 . Va—Vo , Vat Vo — 2V5 


h2 © Qph2 h2 


Zz 


We don’t need to simplify this equation in step 7: 


Vot+Ve—2V5 . Va—Vo | Vat Vo — 2V5 


h2 © Qph2 h2 


r z 


Let us find Vs 


V5 =? 
Let us multiply both sides 


| -2ph2h? 
We obtain 


32 


(Se) _ we 4K -% 
Ps 


(19.1) 


(19.2) 


(19.3) 


(19.4) 


(19.5) 


(19.6) 


(19.7) 


(19.8) 


(19.9) 


430 


432 


433 


434 


435 


436 


437 


438 


439 


442 


QpVeh2 + 2pVgh2 — 4pVsh2 + Veh? — Voh?2+ 
+2pVih2 + 2pVeh2 — 4pVsh2 = 0 


Let us simplify this equation: 


Vs (4ph2 + 4ph2) = V2(2ph2 — h2) + Va(2ph2 + h2)+ 
+2ph2V4 + 2ph2Ve 


So we have: 


V54p(h2 + hz) = Voh2 (2p — 1) + Veh? (2p + 1) + Va2phi, + Ve2ph; 


19.3 Final forms of relaxation formula 


19.3.1. zrLV_RELAX5_P5 A 


hz & be 


hz (2p — 1) hz (2p + 1) hy 
Vo= Go(hE EZ)? dp(h2 eB) © ORE A 


19.3.2 zrLV_RELAX5_ P5_B 


This formula is identica to formula A: 


pe aie 
h2(Qp—1),, | he(Qp+1) h2 


+ acne nay" * acne ny (MM) 


> Ap(h? +h?) 


19.3.3 zrLV_RELAX5_P5_C 


h,=h,py=h 
ops 2p+1 1 
Vs = Vo 4 Vg 4 Va + Ve 
5 8p 2 8p 8 qa 6) 


19.3.4 zrLV_RELAX5 P5 D 


This formula is identical to formula C: 


h,=h,p=h 
Opt 2p+1 1 
Vs = Vo 4 Ve 4 Va + Ve 
5 8p 2 am 8 qiV4+ 6) 


(19.10) 


(19.11) 


(19.12) 


(19.13) 


(19.14) 


(19.15) 


(19.16) 


443 


444 


445 


446 


447 


448 


449 


450 


451 


452 


453 


454 


20 Relaxation formula for node P6 
20.1 Node description 


Right edge of mesh ZR. 


20.2 Calculation of relaxation formula 


Laplace equation at node Ps 


ONG af ONO: oa HO Mei \ 23% 
Or? Et Op? Pes: Oz? 
Po 6 Pe 
Approximation of partial derivatives of V(,,) at node Ps 


(note: r = ph, (such as in Pierre Grivet’s book) ) 


(Se) _ Be 4+ 2M% 
Po 


Or? hr h2 


Ue 


LOVig35 6 VO NG Vo — V3 
Bop? le 2ph?2 


r Or - 


f G6z+ 


2 Vez4+—Ve Ve—Vs 
(° Vier) \ he hz _ Vs — Vo 
Oz? = 7 
Po 


he We chy 
Let us substitute approximations to Laplace equation. 


V3+V9—2Ve _ Vo—V3 Vs—Ve | Ge24 
h2 © Qph2 ° hh, 


r. 


=0 


V3 +Vo—2V6 | Vo— V3) Vs— Ve _ _ 9o24 


; 2ph?. he hz 


Let us find Vg 


Let us multiply both sides 


We obtain 


34 


(20.1) 


(20.2) 


(20.3) 


(20.4) 


(20.5) 


(20.6) 


(20.7) 


(20.8) 


(20.9) 


2pV3hz + 2pVoh2 — ApVeh2 + Voh2 — V3h2+ 


20.10 
+2pVsh2 — 2pVeh2 = —2pgez4+hzh2 ove 


Let us simplify this equation: 


Vo (4ph2 + 2ph2) = V3(2ph2 — h2) + Vo(2ph2 + h2) + Vs2ph2+ 


2pgez+hzh> a 


So we have: 


Vo2p(2h2+h7) = Vah2(2p—1)+Voh2(2p+1)+Vs2ph7+2pge.4hzh7 (20.12) 


20.3 Final forms of relaxation formula 


20.3.1 zrLV_RELAX5 P6_A 


hz F hp 
96z+ F# 0 

h2(2p —1 h2(2p+1 hi ephgh? 

6 = ee ) 34 ee } 94 2 3V5 + Bet 2 (20.13) 
2p(2h2 + h2) 2p(2h2 + h2) 2h2 + h2 2h2 + h2 
20.3.2 zrLV_RELAX5 P6 B 
hz Fh, 
96z+- = 9 
hQp—A h2(2p+1 h? 


~ Qp(2h2 + h2) 2" Qp(2h2 +h2) °°" 2h2 + h2 


20.3.3 zrLV_RELAX5_P6_C 


h,=h,p=h 
96z+ Fx 0 
pete, Sop eeae, dee 2 Hee 
Ve= V3 4 Vo + =V5 4 20.15 
6 ep Gp ot 3% 3 ( ) 
20.3.4 zrLV_RELAX5 P6 D 
ee a 
962+ = 9 
2p —1 2 1 1 
Vee eee (20.16) 
6p 6p 3 


466 


467 


468 


469 


470 


472 


473 


474 


475 


476 


21 Relaxation formula for node P7 
21.1 Node description 


Left, upper corner of mesh ZR. 


21.2 Calculation of relaxation formula 


Laplace equation at node P; 


OPViag 10ViGr PV er 
eee ee 
P; Be P; 
Approximation of partial derivatives of V(,,,) at node P; 


(note: r = ph, (such as in Pierre Grivet’s book) ) 


f Gir+ 


~ 


OView) _ eee oe 2 WA = Vy 
Or? hy hz hp 


(CA) SS, 1 Vor+ bee V4 _ V7 ay Grr+hr a V4 _ 
Pr 


r Or op: Dp 2ph?2 
Vr -—Va | Gtr+ 
Qph2 ph, 
PV) a es an _ VV ote 
Oz? h, he hz 
Pr; 


Let us substitute approximations to Laplace equation. 


Va—Ve | Gtr+ | Va-— Va, Gtr+ | Va -— Ve 97z- 


T t =0 
h2 hy Qph2 * 2h, h2 h. 
Va-Ve  Ve- Va | Va Ve Gtr Stet Ite 
h2 Qph2 h2 hig Dp Pig 


Let us find V7 


Let us multiply both sides 


36 


(21.1) 


(21.2) 


(21.3) 


(21.4) 


(21.5) 


(21.6) 


(21.7) 


(21.8) 


477 


478 


480 


481 


482 


483 


We obtain 


QpVih2 


| -2ph2h? 


— 2pV7h2 + Vrhz — Vah2 + 2pVeh? — 2pVzh? = 
—2pg7r+h2hy = Grr help AF 2pg7z—-hzh? 


Let us simplify this equation: 


Vz (2ph2 — h2 + 2ph?) = Va(Qph2 — h2) + Va(Qph7Z)+ 
2pgrrhzhy + Grr4hzhy — 2pg7z—hzh? 


So we have: 


V7((2p — 1)h2 + 2ph2) = Vih2(2p — 1) + Ve2ph2+ 
2pg7r+hzhy. + Grr hey = 2p972— hzh? 


21.3 Final forms of relaxation formula 


21.3.1. zrLV_RELAX5 P7_A 


V= 


hz F hy 
I7z- # 0 
Gir+ = 0 


h2(2p — 1) 


ais Vet 


(2p — 1)h2 + 2ph2* 


(2p — 1)h2 + 2ph2 


(2p a8 1)g7r+h2hy _ 2p972— hh? 


(2p — 1)h2 + 2ph? 


21.3.2 zrLV_RELAX5 P7_B 


V= 


hz # hy 
97z- = 0 
G7r+ = 90 


h2(2p — 1) 


2ph? 


(2p — 1)h? + 2ph?* 


(2p — 1)h2 + 2ph? Ys 


37 


(21.9) 


(21.10) 


(21.11) 


(21.12) 


(21.13) 


(21.14) 


484 


485 


486 


487 


21.3.3 zrLV_RELAX5 P7_C 


ep al 


V7 


hz =h,=h 
97z—- FO 
Grr+ FO 


2p h((2p + 1) 97+. — 972-) 


~ Ap—1 2 


Vg + 


4p —1 


21.3.4 zrLV_RELAX5 P7_D 


4p — 1 


38 


(21.15) 


(21.16) 


488 


493 


494 


495 


496 


497 


498 


22 Relaxation formula for node P8& 
22.1 Node description 


Upper edge of mesh ZR. 


22.2 Calculation of relaxation formula 


Laplace equation at node Ps 


ON GO af ONO, 5 HO Ven \, 23% 
Or? Gt Oe ss Oz? 
Ps 8 Ps 


Approximation of partial derivatives of V(,,) at node Pg 


(note: r = ph, (such as in Pierre Grivet’s book) ) 


Q 


Or? 


Var —V5 Vg—V5 
(Se) he _ Va Ve | g8r+ 
Ps 


(1%¥en) lL Verg — V5 Va + gsr+hr — V5 _ 
PB: 


hy. hz hp 


r Or op Dp 2ph2 
Vg—V5s | 98r+ 
Qph2 ph, 


Oz? hz h2 


Zz 


(Se) gti ane GN 
Ps 


Let us substitute approximations to Laplace equation. 


V5—Ve  gsrt . V8—V5 . gart . Vr + Vo — 2V2 


=0 
h2 le 2ph?2 2ph, h2 
V5—Va | Vs—Vs | Vit Vo~2V8 _  98r+  98rt 
h2- ' Qh? h2 hie > Dp 
Let us find Vg 
Vg =? 


Let us multiply both sides 


39 


(22.1) 


(22.2) 


(22.3) 


(22.4) 


(22.5) 


(22.6) 


(22.7) 


(22.8) 


| -2ph2h? (22.9) 
499 We obtain 


2pVsh2 — 2pVgh2 + Veh2 — Vsh2 + 2pVehz + 2pVoh; — 4pVgh?, = 


(22.10) 
—2Wggr4 herby — ger+hzhy 
500 Let us simplify this equation: 
Vg (2ph2 — h2 + Aph?) = V5(2phz — h2) + (Vz + Vo)2ph2+ ipod 
2pger+h2hy ae gsr+h2h, , 
501 So we have: 
Va((2p — 1)h2 + 4ph?) = V5h2(2p — 1) + (V7 + Vo)2ph2+ eae 
2pger+h2hy. oy gar+ hh, 
52 22.30 Final forms of relaxation formula 
503 22.3.1 zrLV_RELAX5 P8 A 
hee he 
gsr+ FO 
= h2(2p — 1) 2ph2 
R= z = V. Tr 
° (Qp— Dh? + 4ph2 is (2p — 1)h2 + 4ph?2 vee (22.13) 


(2p + 1h2hy Gert 
(2p — 1)h2 + 4ph?2 


sos 22.3.2 zrLV_RELAX5 P8 B 
hz & Np 
98r+ = 0 
h2(2p — 1) af 2phr 
(2p — 1)h2 + 4ph2° " (2p — 1)h?2 + 4ph? 


506 


Vg= (Vz + Vo) (22.14) 


507 22.3.3 ZrLV_RELAX5 P8_C 


hz=h-=h 
gsr+ FO 
508 op—1 Op 
V3 = V5 4 (V7 + Vo)+ 
ie Sa, (22.15) 
(2p + 1)hgsr+ 
6p —1 


40 


soo 22.3.4 zrLV_RELAX5 P8 D 


510 


41 


(22.16) 


516 


517 


518 


519 


520 


23 Relaxation formula for node P9 
23.1 Node description 


Right, upper corner of mesh ZR. 


23.2 Calculation of relaxation formula 


Laplace equation at node Po 


OPViag 10ViGr PV er 
Gn) (1Mem\) (PV en) 
Or? & r Or Po Oz? 

9 


(note: r = ph, (such as in Pierre Grivet’s book) ) 


vs 
Po 


Approximation of partial derivatives of V(,,-) at node Po 


99r+ 


d 


Or? hy h?2 


2 Vor+—Vo Vo—Vo 
€ se | he hr _ Ve— Vo 
P9 


hy 


(CA) , 1 Vor+ — Vo = Vo + Jor+Rr — V6 = 
Po 


r Or Se Dlg 2ph2 
Vo-—Ve | Gort 
Qph2° 2ph> 
Voz —Vv 97 
P Vz) ae ; i : tte = Joz— | Vg — Vo 
Az? h. lin? ie 
Po 


Let us substitute approximations to Laplace equation. 


Ve—Vo . gor+ , Vo—Vo | Gort , V8—Vo , Gozt _ 


h2 hp  2ph2 © phy R20 hy ° 
Ve—Vo | Vo-Ve | Va-Vo _ ort 99rt 9924 
h2 2ph2 h2 hig Dpiy fig 


Let us find Vo 


Let us multiply both sides 


42 


(23.1) 


(23.2) 


(23.3) 


(23.4) 


(23.5) 


(23.6) 


(23.7) 


(23.8) 


522 


523 


525 


526 


527 


528 


We obtain 


2pVeh2 


| -2ph2h? 


— 2pVah? + Voh2 — Veh2 + 2pVgh? — 2pVoh? = 
—2pgor+h2hy oa gorth2hy = 2pgoz+hzh2 


Let us simplify this equation: 


Vo (2ph2 — h? + 2ph?) = V6(2ph2 — h2) + Vg(Qph2)+ 
2pgor+hzhy a gore heh, or 2pgor+hzh? 


So we have: 


Vo ((2p — 1)h2 + 2ph2) = Veh2 (2p — 1) + Ve2ph?2+ 
(2p a 1) gor+h2h, a 2pgoz+hzh? 


23.3 Final forms of relaxation formula 


23.3.1 zrLV_RELAX5 P9_A 


Y= 


hz # hr 


gJ9 


_ £0 


Jort+ = 0 


h2(2p — 1) 


2ph? 


(2p — 1)h2 + 2ph2 ° 


(2p + 1)gor+h2hy 4 


(2p — 1)h2 + 2ph2 


Vet 


+ 2pgg2+hzh? 


23.3.2 zrLV_RELAX5 P9_B 


Vg = 


(2p — 1)hZ + 2phe 


he Rp 
99z- = 0 
Jor+ = 0 


h2(2p — 1) 


2ph? 


(2p — 1)h? + 2ph?® 


(2p — 1)h2 + 2ph2 


43 


Vg 


(23.9) 


(23.10) 


(23.11) 


(23.12) 


(23.13) 


(23.14) 


529 


530 


532 


533 


551 


552 


553 


23.3.3 zrLV_RELAX5 P9 C 
hz, =h,=h 
Joz— FO 
Jor+ #0 


2p h((2p + 1)gor+ + 992+) 
Ve V; 
" 4p —1 a 4p — 1 


_ 2p—1 
~ Ap—1 


Vo (23.15) 


23.3.4 zrLV_RELAX5 P9 D 


Ve+ Vg (23.16) 


References 


[1] P. Grivet, Electron Optics, Second (revised) English edition. Pergamon 
Press Ltd., 1972. 


[2] J. R. Nagel, “Solving the generalized poisson equation using the fi- 


, _hite - difference method (fdm)..’ https: //my.ece.utah.edu/~ece6340/ 
LECTURES/Feb1/, 2012. [Online; accessed 2-February-2023]. 


[3] kryomaxim, “taylor expansion in cylindrical coordinates.” 


,  jattps://math.stackexchange .com/questions/1133311/ 


taylor-expansion-in-cylindrical-coordinates, 2015. [Online; 
accessed 9-May-2023]. 


[4] A. Septier(ed.), Focusing of Charged Paticles. Volume I. New York and 
London, Academic Press, 1967. 


[5] A. Septier(ed.), Applied Charged Paticle Optics, part A. New York and Lon- 
don, Academic Press, 1980. 


[6] D. W. O. Heddle, Electrostatic Lens Systems. Second Edition. \nstitute of 
Physics Publishing, Bristol and Philadelphia, 2000. 


[7] B. Paszkowski, Optyka Elektronowa, wydanie II, poprawione i uzupetnione. 
Panstwowe Wydawnictwa Naukowo - Techniczne, Warszawa, 1965. 


[8] B. Paszkowski, Electron Optics [by] B.Paszkowski. Translated from the Pol- 
ish by George Lepa. English translation edited by R. C. G. Leckey. London, 
lliffe; New York, American Elsevier Publishing Company Inc., 1968. 


44 


